14  冗余分析

冗余分析（redundancy analysis，RDA）是一种回归分析结合主成分分析的排序方法，也是多因变量（multiresponse）回归分析的拓展。从概念上讲，RDA是因变量矩阵与解释变量之间多元多重线性回归的拟合值矩阵的PCA分析。

14.1 示例1

14.1.1 首先加载R包

library(knitr)
library(tidyverse)
library(agricolae)
library(vegan)
library(ggforce)

14.1.2 加载数据集

以某处试验田数据为例，拟冗余分析明确稻田N_2O排放同土壤NH_4^+-N、土壤NO_3^--N、氧化还原电位Eh、土壤pH、地上部分氮素提取N_{updake}和产量Yield的关系和各因素的贡献。具体数据集见表 14.1。

NH4+-N

dt <- readxl::read_xlsx(here::here("RES/dataset/PCA.xlsx"),sheet="11")
dt.data <- dt %>% select(-trt)
dt %>% kable(digits = 2)
表 14.1 原始数据集
trt N2OEmissions N_NH4 N_NO3 Eh pH Tem Nuptake Yield
ICFZ0 1.31 11.28 1.37 358.48 6.14 23.33 99.77 7.42
ICFZ0 1.00 12.60 1.51 337.21 6.13 23.88 91.09 7.49
ICFZ0 1.06 12.69 1.26 388.40 6.11 23.14 94.94 7.27
IAWDZ0 2.84 12.90 1.48 415.84 6.21 23.40 103.36 7.16
IAWDZ0 3.08 9.90 1.29 426.13 6.19 24.02 102.96 8.39
IAWDZ0 2.92 11.15 1.54 423.48 6.19 22.98 92.43 8.09
IAWDZ10 1.91 17.27 1.83 371.37 6.29 23.33 119.85 9.82
IAWDZ10 2.07 14.68 1.44 376.14 6.26 24.01 120.78 8.89
IAWDZ10 2.14 14.98 1.69 351.42 6.25 23.14 123.49 9.97
ICFZ0 0.54 8.84 2.37 159.19 6.43 20.65 94.33 7.99
ICFZ0 0.41 8.99 2.48 140.97 6.49 21.31 96.10 7.24
ICFZ0 0.54 8.94 2.25 146.49 6.47 20.29 98.09 7.16
IAWDZ0 0.87 8.54 2.51 271.21 6.54 21.42 98.16 7.80
IAWDZ0 0.86 8.14 2.34 216.85 6.57 21.32 107.03 7.61
IAWDZ0 0.79 7.89 2.49 228.48 6.55 20.44 106.04 7.70
IAWDZ10 0.49 10.86 2.55 215.49 6.58 21.41 108.98 8.73
IAWDZ10 0.60 11.51 2.89 199.90 6.58 21.21 112.19 8.90
IAWDZ10 0.61 10.40 2.46 212.04 6.60 20.57 118.03 8.79

14.1.3 数据检测

这样可以得到环境因子数据集为表 14.2。

data.X <- dt.data  %>% select(-N2OEmissions,-Yield)
data.Y <- dt.data %>% select(N2OEmissions,Yield)
data.X %>% kable(digits = 2)
表 14.2 环境因子数据集
N_NH4 N_NO3 Eh pH Tem Nuptake
11.28 1.37 358.48 6.14 23.33 99.77
12.60 1.51 337.21 6.13 23.88 91.09
12.69 1.26 388.40 6.11 23.14 94.94
12.90 1.48 415.84 6.21 23.40 103.36
9.90 1.29 426.13 6.19 24.02 102.96
11.15 1.54 423.48 6.19 22.98 92.43
17.27 1.83 371.37 6.29 23.33 119.85
14.68 1.44 376.14 6.26 24.01 120.78
14.98 1.69 351.42 6.25 23.14 123.49
8.84 2.37 159.19 6.43 20.65 94.33
8.99 2.48 140.97 6.49 21.31 96.10
8.94 2.25 146.49 6.47 20.29 98.09
8.54 2.51 271.21 6.54 21.42 98.16
8.14 2.34 216.85 6.57 21.32 107.03
7.89 2.49 228.48 6.55 20.44 106.04
10.86 2.55 215.49 6.58 21.41 108.98
11.51 2.89 199.90 6.58 21.21 112.19
10.40 2.46 212.04 6.60 20.57 118.03

14.1.3.1 单峰模型或线性模型判定

根据DCA函数的结果，可判定采用CCA分析，还是RDA分析。

decorana(data.Y)

Call:
decorana(veg = data.Y) 

Detrended correspondence analysis with 26 segments.
Rescaling of axes with 4 iterations.
Total inertia (scaled Chi-square): 0.0481 

                        DCA1      DCA2      DCA3      DCA4
Eigenvalues          0.04808 0.0480779 4.808e-02 4.808e-02
Additive Eigenvalues 0.04808 0.0000000 0.000e+00 0.000e+00
Decorana values      0.04808 0.0002719 6.603e-05 3.815e-05
Axis lengths         0.32675 0.3267544 3.268e-01 3.268e-01
判断原则

判断原则

如果DCA1的Axis Lengths>4.0,就应选CCA（基于单峰模型，典范对应分析)；如果在3.0-4.0之间，选RDA和CCA均可；如果小于3.0, RDA的结果会更合理（基于线性模型，冗余分析）

根据DCA1列中Axis Lengths为0.099244,可知Axis Lengths<3.0, 即采用RDA更为合适。

14.1.3.2 冗余分析

result <- rda(data.X ~ N2OEmissions, data=data.Y, scale = TRUE)
result2 <- result %>% summary()

但Explained components只有1项，不能作图。等下次有时间了再找数据。

XX <- envfit(result,data.X)

dtRDA <- data.frame(r=XX$vectors$r) %>% rownames_to_column(var="variable")
Explained <- result2$cont$importance["Proportion Explained",1:2] %>% sum()
UnExpalined <- 1 - Explained

dtExplained <- data.frame(Category = forcats::fct_inorder(c("Explained", "Unexplained")),
                          contrib = c(Explained, UnExpalined))
compo <- XX$vectors$arrows %>% as.data.frame()
r2 <- XX$vectors$r 
p <- XX$vectors$pvals
compo %>% add_column(r2,p)
表 14.3 环境因子对N2O排放的贡献
RDA1 PC1 r2 p
N_NH4 0.9933848 -0.1148333 0.5420960 0.003
N_NO3 -0.9908248 0.1351526 0.8976482 0.001
Eh 0.9630207 0.2694275 0.9364041 0.001
pH -0.9298926 0.3678311 0.9012294 0.001
Tem 0.9997646 -0.0216944 0.9372436 0.001
Nuptake 0.5145326 0.8574708 0.1609866 0.267 
dtRDA$cumSum <- cumsum(dtRDA$r/sum(dtRDA$r)*Explained)
dtRDA$Category <- "Explained"
ggplot()+
  geom_bar(data=dtExplained,aes(x=contrib,y="",fill=Category),stat="identity")+
  geom_text(data=dtRDA,aes(x=cumSum+UnExpalined-r/sum(r)/2,1.5+0.3,label=sprintf("%.1f%%",r*100)),check_overlap=TRUE)+
  geom_text(data=dtRDA,aes(x=cumSum+UnExpalined-r/sum(r)/2,1.5+0.6,label=variable),vjust = -0.2,check_overlap=TRUE)+
  geom_ellipse(data = dtRDA,aes(x0=cumSum+UnExpalined-r/sum(r)/2,y0=1.5+0.6,a=0.05,b=0.1,angle=0))+
  geom_segment(
    data = dtRDA,
    aes(
      x = cumSum+UnExpalined-r/sum(r)/2,
      y = 1.5+0.6,
      xend = cumSum +UnExpalined-r/sum(r)/2,
      yend =1.5
    ),
    linewidth = 0.25,
    arrow = arrow(
      angle = 15,
      length = unit(0.25, "cm"),
      type = "closed"
    ),
    linetype = 1,
    # size = 0.6,
    colour = "red"
  )+
  coord_polar(theta="x",direction = -1,start = Explained/2*6.4)+
  theme_void()
图 14.1 环境因子对N2O方差变化贡献